r"""Solver for linear systems of the form:

.. math::

    \frac{\partial Q}{\partial t} + A\frac{\partial Q}{\partial x} = 0

The Jacobian matrix :math:`A` is a constant.

"""

# base solver module
import solver
import flux.acoustics.acoustics_flux as flux

import numpy
from os.path import join

class LinearSolver(solver.Solver1D):
    r"""Solver for the linear systems.

    **The numerical Method**

    The numerical update formula used is the fluctuation form for the
    cell updates as described by LeVeque (section 6.13 ):

    .. math::

        Q_{i}^{n+1} = Q_{i} - \frac{\Delta t}{\Delta x} \left(
        A^{+}\Delta Q_{i-\frac{1}{2}} + A^{-}\Delta Q_{i +
        \frac{1}{2}} \right) - \frac{\Delta t}{\Delta
        x}\left(\tilde{F}_{i + \frac{1}{2}} - \tilde{F}_{i -
        \frac{1}{2}} \right ),

    where,

    .. math::
        
        A^{+}\Delta Q_{i - \frac{1}{2}} = \sum_{p=1}^{m}
        \lambda^{+}W_{i - \frac{1}{2}}^{p}

        A^{-}\Delta Q_{i - \frac{1}{2}} = \sum_{p=1}^{m}
        \lambda^{-}W_{i - \frac{1}{2}}^{p}

        \tilde{F}_{i - \frac{1}{2}} = \frac{1}{2}\sum_{p=1}^{m}
        |\lambda^{p}| \, \left( 1 - \frac{\Delta t}{\Delta
        x}|\lambda^{p}|\right) \tilde{W}_{i-\frac{1}{2}}^{p}

    The first term in brackets in the update formula corresponds to
    the standard Godunov's upwind flux. The second term provides the
    high resolution corrections to make the resulting scheme formally
    second order. To avoid numerical oscillations that result from
    second order schemes, TVD limiters are used to limit the waves in
    the Riemann solution:

    .. math::

        \tilde{W}_{i-\frac{1}{2}}^p =
        \alpha_{i-\frac{1}{2}}^{p}\,
        \phi(\theta_{i-\frac{1}{2}}^{p}),

    where the limiter function is given by

    .. math::

        \theta_{i - \frac{1}{2}}^{p} = \frac{
        \alpha_{I-\frac{1}{2}}^{p} }{ \alpha_{i-\frac{1}{2}}^{p} }


        with\,\,  I = \begin{cases}
              i - 1 & if \lambda^p > 0 \\
              i + 1 & if \lambda^p < 0
              \end{cases}

    This flux limiter method has a five point stencil as opposed to
    the three point stencil used by the Lax Wendroff method.

    Reference reading for schemes of this form is found in LeVeque.

    **Implementation**

    We choose to implement all schemes in PyHCL in the standard
    conservative form:

    .. math::

        Q_{i}^{n+1} = Q_{i}^{n} - \frac{ \Delta t }{ \Delta x }\left[
        F_{i+\frac{1}{2}} - F_{i - \frac{1}{2}} \right ]


    We compute the flux for each face of a numerical cell. In the 1D
    case, there are only two faces at :math:`i + \frac{1}{2}` and
    :math:`i - \frac{1}{2}`.Individual flux functions are responsible
    for determining the flux at the cell interfaces. Refer to
    :file:`solvers.flux.acoustics.acoustics_flux.py` for an example of
    a linear flux function.
    
    """
    def __init__(self, nvar, tf, grid=None, flux_type=0):
        solver.Solver1D.__init__(self, nvar, tf, grid=grid)
        self.defaults = dict( flux_type=flux_type ) 
    
    def compute_dt(self):
        r"""Compute the time step.

        For a uniform grid of size :math:`\Delta x` and a CFL number
        :math:`c`, the time step is given by :math:`\Delta t = c
        \Delta x`

        """
        return self.cfl * self.grid.dx
